<html>
    <head>
        <meta charset="UTF-8">
		<style>
			body{
				background:#ecfbdb url('../res/bg.bmp')
			}
			#wmd-preview{
				margin: 6% 13%;
				padding: 9% 3%;
				background:#fcf9f3;
			}
			#program-window{
				background-color:#679;
				padding:40px;
				float:right;
				z-index: 9999;
				position: relative;
			}
		</style>
        <script type="text/javascript" src="../lib/mathjs/9.4.4/math.min.js"></script>
		<script type="text/javascript" src="./math.js"></script>
		<script type="text/javascript" src="./mathTest.js"></script>
		
        <script>
			var cvx = null, ctx = null, pts = [];
			var pts = [];

			function exec() {
				//画点
				ctx.clearRect(0, 0, cvs.width, cvs.height);
				pts.forEach(function(value) {
					ctx.fillStyle = "green";
					ctx.beginPath();
					ctx.arc(value.x, value.y, 3, 0, 6.3);
					ctx.fill();
				});

				if (pts.length <= 1)
					return;

				let x_y_1 = math.zeros(pts.length, 3)
				for(let i = 0; i < pts.length; i++)
					x_y_1.subset(math.index(i, [0, 1, 2]), [pts[i].x, pts[i].y, 1])
				res = math.eigens(math.multiply(math.transpose(x_y_1), x_y_1))
				res = res[1]
				a = res.get([2, 0])
				b = res.get([2, 1])
				c = res.get([2, 2])
				
				ctx.beginPath()
				ctx.strokeStyle = 'orange'
				ctx.moveTo(-a * c / (a**2 + b**2) + 100000 * (-b), -b * c / (a**2 + b**2) + 100000* a)
				ctx.lineTo(-a * c / (a**2 + b**2) - 100000 * (-b), -b * c / (a**2 + b**2) - 100000 * a)
				ctx.stroke()
			}

			function getEventPosition(ev) {
				var x, y;
				if (ev.layerX || ev.layerX == 0) {
					x = ev.layerX;
					y = ev.layerY;
				} else if (ev.offsetX || ev.offsetX == 0) {
					// Opera
					x = ev.offsetX;
					y = ev.offsetY;
				}
				return {
					x: x,
					y: y
				};
			}

			function on_btn_clear_clicked() {
				pts = [];
			}
			
			function onload() {
				cvs = document.getElementById('canvas');
				ctx = cvs.getContext('2d');
				cvs.addEventListener('click', function(e) {
					let pos = getEventPosition(e);
					pts.push(pos);
				}, false);
				window.setInterval("exec()", 100);
				console.log("程序设置完成！")
			}
        </script>
		<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
		<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
    </head>
    <body onload="onload()">
		<div id="wmd-preview" class="wmd-preview">
			<div id="program-window">
				<button onclick="on_btn_clear_clicked()">清空</button>
				<canvas style="border:2px solid #eee; background-color:white" width="314" height="200" id="canvas">浏览器不兼容。</canvas>
				<div>
					<br>
					Powered by:<br>
					<li><a href="http://mathjs.org/index.html">Math.js</a>-<a href="http://www.apache.org/licenses/LICENSE-2.0">Apache 2.0 License</a></li>
				</div>
			</div>
			<p data-anchor-id="vsa5">
			\(\newcommand\v[1]{\boldsymbol{#1}}\)
			<p>
			笔者设计了一个拟合直线的程序，放在文章的右边。阅读本文前，读者不妨试着点击鼠标，在图中添加几个点，提前获得关于最小二乘法的感性认识。
			</p>
			<p>
			直线的拟合问题是运用最小二乘法的一个很好的范例，我们就利用这个问题来理解最小二乘法。</p>
			<h2>问题</h2>
			已知二维平面上的一些列点\((x_k, y_k)\)，\(k=1,2,3,\cdots\)要求用最小二乘法拟合直线。</p><p>
			<h2>解</h2>
			设直线方程为\(ax+by+c=0\)，则有$$\left(\begin{matrix}
			x_1 & y_1 & 1\\ 
			x_2 & y_2 & 1\\
			\vdots & \vdots & \vdots\end{matrix}\right)
			\left(\begin{matrix}
			a\\b\\c
			\end{matrix}\right)=0$$<br>
			方程满足\(A\v{x}=0\)的形式。其中\(A\)是\(k\times 3\)的矩阵，\(x=(a, b, c)^T\). <br>
			所谓最小二乘法，即是求解\(x\)使得$$||Ax||^2=(Ax)^T Ax=x^T A^T A x$$的值最小。<br>
			这个问题显然有一个解\(x=0\)，但这个解没有意义，因为它不能构成直线的系数，人们把这个解叫做“平凡解（trivial solution）”。
			<br>
			为了排除平凡解，引入一个约束条件：
			$$x^Tx=1$$
			<br>
			问题变为约束条件下求最小值的问题，可以运用拉格朗日乘数法。令 
			$$L(x)=x^TA^TAx-\lambda(x^Tx-1)$$ <br>
			\(L(x)\)最小时，\({\partial L(x)\over \partial x} = 0\)，因此
			$$A^TAx - \lambda x = 0$$ <br>
			因此\(\lambda\)是\(A^TA\)的特征值，\(x\)是\(A^T A\)对应\(\lambda\)的一个单位特征向量，记\(x=\boldsymbol{e_\lambda}\)。代入得\(L(\boldsymbol{e_\lambda})=\lambda \boldsymbol{e_\lambda}^T \boldsymbol{e_\lambda}=\lambda\)。\(\lambda\)最小时，\(L(x)\)取得最小值。换言之，\(A^TA\)的最小特征值对应的单位特征向量即是所求的解。</p>
			<p>
			求矩阵的特征特征向量有几种方法。注意这里的\(A^T A\)是实对称矩阵，可以用<a href="https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm">Jacobi eigenvalue algorithm</a>求特征值。更常见的方法则是使用SVD分解。
			</p>
			<h3>SVD分解求\(A^T A\)的特征值</h3>
			根据SVD分解，
			$$A=UDV^T=
			\left[\begin{matrix}
			\boldsymbol{u_0} | \cdots | \boldsymbol{u_{p-1}}
			\end{matrix}\right]
			\left[\begin{matrix}
			\boldsymbol{\sigma_0} & & \\
			 & \ddots & \\
			  & & \boldsymbol{\sigma_{p-1}}
			  \end{matrix}\right]
			\left[\begin{matrix}
			\boldsymbol{v_0^T}\\ \hline
			\vdots\\ \hline
			\boldsymbol{v_{p-1}^T}
			\end{matrix}\right]
			$$
			注意到\(A^TA=(UDV^T)^T(UDV^T)=VDU^TUDV^T=V\Lambda V^T\)
			其中，\(\Lambda\)的第i行的非0元素\(\lambda_i=\sigma_i^2\)
			两边右乘\(V\)得到
			$$(A^TA)V=V\Lambda$$
			即是说，\(V\)的每一个列向量\(\boldsymbol{v_j}\)都是\(A^TA\)的特征向量，对应的特征值是\(\lambda_j\).
			<h2>总结</h2>
			<p>為了解決問題，筆者引入約束條件\(x^Tx=1\)，然後運用拉格朗日乘數法來求解最小值。實際上，還有其它種形式的約束條件，例如，可以約定\(b=1\)，但這樣一來，
			直線方程從一般式退化成斜截式，就無法表示平行于y軸的直線了。儘管如此，在工程中人們還是常常運用斜截式。
			</p>
			<p>
			若是採用斜截式，經過求導等步驟，一樣可以求出直線的參數。但不論如何，最小二乘法是使殘差平方和最小化的方法，換言之，只要能使殘差平方和最小化，無論使用何種方法來求取x的值，都可以叫做最小二乘法。</p>
			<p>
			最小二乘法不僅可以用於直線的擬合，還可以用於擬合圓、二項式、平面……等等模型。歷史上，高斯就成功用最小二乘法預測了谷神星的軌跡。
			</p>
			<h2>参考资料</h2>
			<p>
			<li><a href="http://cs-courses.mines.edu/csci508/schedule/03/SVD.pdf">SVD 分解</a></li></p>
		</div>
    </body>
</html>
